Recycling  Propellants  in  Nonpolluting  Supercritical  Fluids: 

Novel  Computational  Chemistry  Models  For  Predicting  Effective  Solvents 

PP-695/E316 

FINAL  TECHNICAL  REPORT 

Betsy  M.  Rice  and  Cary  F.  Chabalowski 
U.  S.  Army  Research  Laboratory 
AMSRL-WM-BD 

Aberdeen  Proving  Ground,  Maryland  21005-5066 


PROJECT  BACKGROUND 

The  conventional  methods  of  demilitarization  of  energetic  components  of  the 
demilitarization  inventory,  open  burning  or  open  detonation  (OB/OD),  are  being  restricted  due  to 
the  resulting  generation  of  pollution.  Additionally,  OB/OD  destroys  valuable  energetic  materials 
that  could  otherwise  be  reused  if  extracted  from  the  stockpiles.  An  alternative  demilitarization 
procedure  to  OB/OD  is  supercritical  fluid  (SF)  extraction  of  the  energetic  materials.  In  addition 
to  eliminating  pollution,  SF  extraction  allows  for  the  recycling  of  the  recovered  ingredients, 
provided  the  process  is  nondestructive.  A  particularly  attractive  nondestructive  solvent  is  carbon 
dioxide,  a  low-cost  environmentally  benign  solvent  with  easily  accessed  critical  parameters. 

Unfortunately,  the  ingredients  in  composite  propellants  and  explosives  have  varying 
degrees  of  solubility  in  SF  CO2,  but  it  has  been  observed  that  solubility  is  enhanced  when  trace 
amounts  of  simple  polar  modifiers  are  added  to  the  SF  solvent.  Morris  et  al.  [1,  2]  reported 
experimental  results  on  the  extraction  of  RDX  from  grains  of  a  nitramine-based  gun  propellant, 
M43,  using  SF  CO2  solvents,  each  of  which  was  modified  by  one  of  34  polar  modifiers.  This 
survey  study  was  performed  in  an  attempt  to  establish  correlations  between  modifier  property  and 
increase  in  the  amount  of  RDX  extracted.  M43  is  composed  of  76%  RDX,  a  16%  mix  of 
polymers  nitrocellulose  and  cellulose  acetate/butyrate  and  8%  plasticizer.  A  grain  of  M43  has  a 
cylindrical  shape  and  has  a  mass  in  the  range  of  1.2- 1.4  g.  The  grains  used  in  the  studies  had  19 
perforations  running  axially,  and  were  coated  with  graphitic  glaze.  In  the  first  study  [1],  at  the 
beginning  of  each  extraction  (41  MPa,  323.2  K),  a  single  grain,  the  CO2  solvent  and  co-solvent 
were  added  to  the  extraction  vessel,  and  the  system  allowed  to  equilibrate  for  two  hours.  The 
extraction  vessel  was  then  flushed  with  neat  CO2.  The  displaced  modified  CO2  solvent  was  then 
expanded  through  a  heated  flow  restrictor  into  an  acetonitrile  solvent  trap,  where  any  extracted 
RDX  was  collected.  Increases  in  the  amount  of  RDX  extracted  in  solvents  containing  polar 
modifiers  relative  to  that  obtained  using  neat  SF  CO2  were  reported  as  extraction  enhancement 
factors  (EEF).  The  results  showed  that  the  most  effective  modifier  is  dimethyl  sulfoxide 
(DMSO),  with  an  EEF  of  69.  Since  the  removal  of  RDX  was  performed  through  a  combined  2- 
hour  static  equilibration  process  followed  by  a  dynamic  extraction,  it  is  possible  that  the 
solubility  limit  had  not  yet  been  reached.  The  solubility  limit  is  defined  as  the  mass  of  a 
particular  solute  dissolved  in  a  mass  of  solvent  that  has  reached  thermodynamic  equilibrium  [1], 
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In  a  case  in  which  the  solubility  limit  has  not  been  reached,  the  EEF  due  to  the  addition  of  polar 
modifier  is  assumed  to  be  a  relative  representation  of  the  enhancement  of  solubility  of  a  solute 
due  to  the  modified  solvent.  Under  such  an  assumption  the  EEF  is  probably  lower  than  the 
enhancement  in  solubility  of  a  solution  in  thermodynamic  equilibrium.  Also,  since  M43 
contains  multiple  components,  it  is  possible  that  the  amount  of  RDX  extracted  is  less  than  that 
which  would  be  extracted  from  a  sample  of  pure  RDX.  In  a  follow-up  study  [2],  Morris  et  al. 
performed  dynamic  SF  extraction  of  RDX  from  M43  propellant  grains  that  had  been  ground.  In 
these  measurements,  Morris  et  al.  report  that  the  extractions  are  still  probably  mass-transfer 
limited,  thus  the  extraction  enhancement  ratios  most  likely  underestimate  the  true  solubility 
enhancement  factors.  However,  Morris  et  al.  report  that  these  measurements  are  closer  to  those 
that  would  correspond  to  the  solubility  limit  than  those  measured  using  bulk  propellant  [1],  The 
amount  of  RDX  extracted  in  this  experiment  was  greater  than  that  extracted  in  the  study  reported 
in  Ref.  1  by  factors  of  ~  8  for  neat  CO2  and  ~4  for  DMSO  and  acetonitrile  modifiers. 

Although  Morris  et  al.  [1,  2]  found  weak  correlations  between  modifier  properties  and 
increase  in  EEF  due  to  the  addition  of  polar  modifiers,  the  dependencies  of  the  enhancement 
factors  on  the  parameters  of  interest  are  not  established.  The  objective  of  this  project  is  to 
determine,  through  the  use  of  well-established  computational  chemistry  techniques,  the  optimal 
physical  conditions  and  chemical  makeup  of  an  effective  SF  CO2  solvent  with  added  polar 
modifier  for  the  extraction  of  a  nitramine  explosive.  This  is  accomplished  with  the  development 
of  accurate  molecular  models  for  use  in  computer  simulation  to  determine  the  conditions  that 
will  result  in  the  maximum  extraction  of  reusable  ingredients  from  waste  energetic  materials.  The 
software  developed  under  this  project  is  designed  to  allow  the  non-expert  user  the  ability  to 
easily  set  up  and  execute  “computer  experiments”  to  determine  the  conditions  associated  with 
optimal  extraction  for  various  chemical  systems.  Additionally,  the  software  can  be  executed  on 
either  Unix-based  or  Windows-based  computer  platforms.  We  have  parameterized 
mathematical  models  for  use  in  the  simulations  to  describe  the  34  co-solvent  systems  measured 
in  the  experimental  study  [1];  these  are  included  in  the  software  available  for  distribution.  We 
also  performed  extensive  analyses  to  understand  the  enhancement  in  the  solubility  of  RDX  due  to 
the  presence  of  three  different  polar  modifiers  in  the  SF  CO2  solvent.  These  studies  were 
centered  on  three  representative  modifiers  whose  EEFs  spanned  the  range  of  the  34  values 
measured  in  the  experimental  study  [1],  We  studied  a  modifier  with  a  large  EEF,  DMSO;  one 
with  a  moderate  EEF,  CH3CN;  and  one  that  had  a  very  small  EEF,  CH3OH.  We  analyzed  the 
dependence  on  solubility  enhancement  due  to  several  parameters,  such  as  concentration,  dipole 
moment,  temperature  and  pressure. 

The  following  describes  the  theory  upon  which  the  determination  of  the  solubility  of  the 
nitramine  explosive  in  the  solvent  is  based,  followed  by  a  description  of  the  models  used  for  the 
determination,  and  a  discussion  of  computational  methodologies.  Results  describing  simulations 
in  which  the  potential  energy  models  were  assessed  and  in  which  solubilities  were  determined  are 
next  described. 

Theoretical  Formulation  for  Calculation  of  the  Solubility 
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When  a  solute  in  the  solid  phase  is  at  equilibrium  with  a  fluid  containing  the  solvent  and 
the  solute  at  temperature  T  and  pressure  P,  the  chemical  potential  of  the  solute  molecule  in  the 
solid  phase,  p  ,  equals  that  in  the  gas  (fluid)  phase,  p  : 

(is(T,P)=  n°(T,P).  (1) 

If  p*  is  the  vapor  pressure  of  the  solute  at  temperature  T,  and  if  p*  is  so  small  that  the 
vapor  at  T  and  p*  can  be  assumed  to  be  an  ideal  gas,  then  one  can  write 

PS(T,  p*)  —  pideaI  gas(T,  p*).  (2) 

For  the  case  when  the  density  of  the  solid  phase  is  assumed  to  be  constant  with  the  variation  of 
pressure  from  p*  to  P  at  a  fixed  temperature,  one  can  write 

ps(T,  P)  =  ps(T,  p*)  +  j  v  dP  =ps(T,  p*)  +  v  (P-p*) ,  (3) 

p * 

where  v  is  the  ratio  of  the  molar  volume  of  the  solid  and  Avogadro’s  number.  The  chemical 
potential  pG(T,  P)  can  then  be  expressed  as  a  sum  of  two  terms:  chemical  potential,  pldeal  gas, 
corresponding  to  the  ideal  gas  at  the  same  temperature  and  pressure,  and  the  residual  chemical 
potential,  pr , 

P°(T,  P)  =  p,deaJ  gas  (T,  P)  +  pr(T,  P)  .  (4) 

Using  Eqs.(l)-(4)  and  the  expression  [3]  for  pldeal  8as  one  can  obtain 

N J  V=  (p*/kT)  exp(-pr/kT)  exp[v(P-p*)/kT] ,  (5) 

where  k  is  the  Boltzmann  constant  and  Ns  is  the  number  of  solute  molecules  in  the  solution 
having  volume  V. 

By  considering  Nc  to  be  the  number  of  CO2  molecules  in  solution  in  volume  V,  the 
solubility  y  can  be  written  as 

y=Ns/Nc=(p*/pkT).exp(-pr/kT)exp[v(P-p*)/kT] ,  (6) 

where  p=  Nc/  V,  the  number  density  of  CO2  molecules.  Equation  (6)  can  be  used  to  compute  the 

solubility  using  the  residual  chemical  potential,  pr.  The  residual  chemical  potential  can  be 
determined  using  Widom’s  test  particle  method  [4]  in  the  NPT  ensemble  for  an  infinitely  dilute 
solution  according  to  the  following  formula  [5] 

pr  -  -kT  In  [  <Vexp(-TVkT)>/<V>].  (7) 

Here,  V  is  the  volume  of  the  fluid,  and  <  >  denotes  the  average  over  all  configurations  of  the 
fluid  and  all  random  orientations  of  the  test  particle.  'F  in  the  above  equation  gives  the  total 
potential  energy  between  a  single  solute  ‘test’  particle  (in  the  absence  of  any  other  solute 
molecules)  and  all  of  the  molecules  in  the  solvent  for  a  given  configuration.  It  may  be  noted  that 
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if  density  p  in  Eq.  (6)  represents  the  number  density  of  solvent  (CC>2+modifier)  molecules  then  y 
would  give  the  ratio  of  number  of  solute  molecules  to  the  number  of  solvent  molecules  [6,  7]. 

Computational  details  of  the  molecular  simulations 

(a)  NPT-MC  simulation  of  solubility 

The  simulation  cell  described  herein  is  composed  of  n  CO2  molecules  and  N  modifier 
molecules  contained  in  a  cubical  box.  Periodic  boundary  conditions  are  imposed  throughout  the 
simulation.  The  values  of  n  and  N  have  been  chosen  according  to  the  desired  molar  concentration 
[N/(n+N)]  of  the  polar  modifier.  The  total  number  of  solvent  molecules  (n+N)  value  are  ~300, 
with  the  exact  value  dependent  on  the  molar  concentration  of  polar  modifier. 

We  have  performed  Monte  Carlo  simulations  (MC)  in  the  isothermal-isobaric  (NPT) 
ensemble.  The  CO2  molecules  are  treated  as  rigid  non-spherical  bodies,  and  the  method  of 
quaternions  [8]  is  used  to  represent  the  orientation  of  the  CO2  molecules  about  the  centers  of 
mass.  One  cycle  of  MC  steps  consists  of  (2n+N+l)  moves:  a  random  move  in  the  orientation  of 
each  of  the  n  CO2  molecules,  (N+n)  random  translational  moves  of  the  centers  of  mass  of  all 
molecules  and  one  random  change  in  the  size  of  the  box.  For  obtaining  the  residual  chemical 
potential  at  one  value  of  T  and  P,  50000-200000  cycles  of  MC  steps  have  been  performed. 

A  50%  acceptance  criterion  is  used  to  select  the  maximum  step  size.  The  steps  sizes  for 
the  box  move,  the  quaternion  move,  and  those  of  the  center  of  mass  moves,  are  found  to  be  close 
to  0.30  A,  0.25,  and  0.40  A  ,  respectively.  The  cutoff  distances  for  the  interaction  potentials  are 
equal  to  one-half  the  size  of  the  simulation  box.  Also,  the  corrections  for  the  long-range 
interaction  potential  have  been  included  [8]. 

(b)  Computation  of  the  residual  chemical  potential 

Equation  (7)  has  been  used  to  compute  the  residual  chemical  potential,  pr .  It  requires  the 
computation  of 'P  and  V  for  a  large  number  of  random  insertions  of  the  test  particle  over  a  large 
number  of  configurations  of  the  CO2  +  modifier  fluid. 

In  order  to  compute  the  potential  'P  for  a  given  configuration  of  the  (n+N)  molecules  in 
the  box  of  volume  V  we  first  select  randomly  the  location  of  the  test  particle  (RDX)  in  the  box. 
We  also  select  randomly  the  orientation  of  the  point  dipole  moment  associated  with  the  test 
particle.  is  then  calculated  using  the  following  equation: 

£  'f'kCCOz)  +2  'Pij(site)  +£  %(dd)  +  'F(long).  (8) 

Here  T'kfCC^)  is  the  interaction  potential  between  the  test  particle  and  k  th  CO2  molecule, 
T'ijfsite)  denotes  the  interaction  energy  between  the  test  particle  and  the  i  th  site  of  the  j  th 
modifier  molecule,  and  T'j(dd)  gives  the  dipole-dipole  interaction  between  j  th  modifier  molecule 
and  the  test  particle  [9].  The  summation  in  this  equation  extends  over  all  molecules  and  groups 
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within  the  cutoff  distance  Rc  from  the  test  particle.  SF(long)  represents  the  usual  long  range 
contribution  of  the  molecules  beyond  the  cutoff  distance.  The  steps  used  to  compute  'F  are 
described  in  greater  detail  in  Ref.  1 1 . 

C.  Calculation  of  solubility  and  estimated  error 

The  solubility  of  RDX  in  polar-modified  CO2  at  supercritical  fluid  temperatures  and 
pressures  has  been  computed  by  using  Eq.  (6)  and  the  residual  chemical  potential  calculated 
using  NPT-MC  methods.  The  partial  density  p  of  CO2  that  is  used  in  Eq.  (6)  was  determined  by 
averaging  the  density  over  all  configurations  of  fluid  throughout  the  MC  simulations.  The  value 
of  v  was  determined  from  the  experimental  density  (1 .806  g/cc  at  300  K)  [12]  and  the  volume 
expansion  coefficient  (0. 191*1  O'3  K'1)  [13]  of  RDX.  The  values  of  p*  have  been  taken  from  the 
literature  [14,  15].  The  enhancement  factor  in  the  solubility  is  defined  as  the  ratio  of  the 
computed  solubility  in  the  polar-modifier  CO2  solvent  to  the  computed  solubility  in  unmodified 
CO2[10], 

The  errors  due  to  the  finite  numbers  of  Monte  Carlo  moves  and  insertions  of  the  test 
particle  procedure  were  estimated  using  the  procedure  described  in  Ref.  10.  This  procedure 
randomly  partitions  the  total  data  generated  in  a  simulation  into  two  sets,  from  which  averages 
are  calculated  and  compared.  This  procedure  is  repeated  hundreds  of  times  for  each  simulation. 
The  estimated  error  factor  in  the  solubility  is  then  determined  in  accordance  with  the  range  of 
variation  in  the  residual  chemical  potential  for  80%  of  the  hundreds  of  sets  so  generated.  For  the 
solubility  results  reported  in  this  paper  the  error  factors  so  estimated  are  in  the  range  1.35  to  1 .95. 
The  computed  results  also  indicate  that  the  error  due  to  the  limited  number  of  insertions  of  the 
test  particle  are  smaller  than  the  error  due  to  the  finite  number  of  MC  moves  for  all  cases.  It 
seems  reasonable  as  the  error  computed  here  due  to  the  finite  number  of  MC  moves  also  includes 
the  error  due  to  the  finite  number  of  insertions. 

Development  of  the  Models 

Analytic  functions  that  describe  interaction  potential  energies  for  chemical  systems  are 
parameterized  using  quantum  mechanical  (QM)  calculations  and  experimental  measurements. 
Assessments  of  the  quality  of  the  potential  energy  function  is  obtained  through  its  subjection  to 
classical  molecular  simulations  [molecular  dynamics  (MD)  or  Monte  Carlo  (MC)].  A  QM 
calculation  is  based  on  fundamental  first  principles,  and  solution  of  the  QM  equations  provides  a 
comprehensive  and  precise  description  of  a  system  of  atoms.  The  information  generated  through 
these  calculations  is  used  to  develop  mathematical  models  of  the  chemical  system;  in  this  case,  the 
solvent  consisting  of  CO2  molecules  and  polar  modifier  molecules,  and  the  solute  molecule.  Once 
these  mathematical  models  are  developed,  they  are  then  used  in  classical  molecular  simulations, 
from  which  time-dependent  or  thermally  averaged  properties  of  the  material  are  determined.  These 
simulations  allow  for  the  performance  of  “computer  experiments”  in  which  the  user  can  specify  the 
temperature  and  pressure,  and  percent  concentration  of  polar  modifier  in  the  solvent.  Although  not 
first  principles,  classical  molecular  simulations  using  realistic  models  result  in  very  accurate 
predictions  of  reaction  dynamics  and  bulk  properties  of  materials.  These  calculations  are  used  to 
determine  the  solubility  of  an  explosive  in  SF  CO2  with  and  without  polar  modifiers. 
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Quantum  Mechanical  Determinations  of  Intermolecular  Interactions 

A  significant  number  of  quantum  mechanical  calculations  were  performed  to  determine  the 
structural  parameters  associated  with  the  equilibrium  conformation  of  the  molecules  in  the  system 
and  interacting  forces  between  the  various  pairs  of  molecules  [16-19].  The  compilation  of  this  data 
allowed  for  the  development  of  the  mathematical  models  that  were  used  in  the  molecular 
simulations,  discussed  hereafter. 

Symmetry-Adapted  Perturbation  Theory  (SAPT)  was  implemented  to  determine  interaction 
energies  between  pairs  of  solute-solvent  molecules  [16-18].  The  focus  of  our  project  is  on  the 
nitramine  explosive  RDX;  however,  it  is  too  large  to  be  effectively  treated  using  the  SAPT 
methodology  and  thus,  a  smaller  representative  chemical  system  must  be  used  in  order  to  determine 
the  interactions  due  to  the  presence  of  a  nitramine  explosive.  For  this  purpose,  we  have  selected  the 
small  molecule,  dimethylnitramine  (DMNA).  RDX  can  be  thought  of  as  a  “trimer”  formed  by  three 
DMNA  molecules  arranged  in  the  ring;  thus,  it  is  expected  that  DMNA  will  have  all  the  required 
characteristic  features  of  the  RDX  molecule  in  terms  of  interactions  with  the  solvent  and  cosolvent 
molecules.  Interactions  between  DMNA  with  CO2,  CH3CN,  CH3OH  and  with  another  DMNA 
molecule  have  been  calculated,  and  analytic  functions  developed  to  describe  these  interactions  as 
functions  of  intemuclear  separation  and  relative  molecular  orientation  [18].  Additionally,  similar 
calculations  have  been  performed  to  develop  CO2-CO2  [17]  and  CH3CN-CO2  [16]  interaction 
potentials  for  use  in  molecular  simulations.  For  complexes  not  involving  DMNA,  the  highest 
available  level  of  S APT  has  been  applied,  while  the  remaining  systems  were  treated  in  a  more 
approximate  manner  due  to  computational  limitations.  These  calculations  indicate  that  favored 
interactions  between  dimer  pairs  are  due  to  an  interplay  between  electrostatic,  dispersion,  induction 
and  exchange  interactions,  none  of  which  can  be  neglected  when  considering  the  shape  and  stability 
of  a  complex.  It  was  found  that  simple  electrostatic  arguments  for  these  systems  lead  to  incorrect 
conclusions  about  the  stabilities  of  complexes,  and  that  often  the  dispersion  contribution  of  the 
interactions  between  the  dimer  pairs  is  as  significant  as  the  electrostatic  interaction.  In  general, 
induction  effects  are  smaller  than  the  other  interactions,  but  in  a  few  cases,  they  are  as  large  as  the 
dispersion  contributions  to  the  energy.  If  the  observations  from  this  work  extend  to  other  systems  of 
large  molecules  (such  as  biological  molecules),  then  the  conventional  interpretations  of  interactions 
based  on  the  electrostatics  alone  cannot  be  trusted  as  reliable. 

A  second  series  of  lower-level  quantum  mechanical  calculations  were  performed  for 
isolated  RDX  molecules  in  various  conformations,  in  order  to  establish  the  electrostatic  potential 
surrounding  the  molecules  [19,  20].  This  information  was  used  to  determine  atom-centered 
partial  charges  that  could  be  used  in  a  standard  Coulombic  model  of  electrostatic  interactions  for 
RDX  with  some  other  molecule.  The  assignment  of  the  electrostatic  charges  was  made  by  using 
the  set  of  atom-centered  monopole  charges  for  the  isolated  molecule  (with  the  structure  fixed  at 
the  experimental  crystallographic  configuration)  that  best  reproduces  the  quantum  mechanically 
derived  electrostatic  potential.  A  variety  of  molecular  properties  were  evaluated  using  the 
quantum  mechanical  information,  including  vibrational  spectra  and  conformational  structures 
and  energies  in  order  to  assess  the  quality  of  the  calculations  used  in  the  model  development 
[19]. 
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Interaction  Potential  Functions  for  use  in  Molecular  Simulation 


RDX/CQ?  systems 


As  indicated  in  the  earlier  section,  computational  limitations  precluded  the  quantum 
mechanical  description  of  interactions  between  pairs  of  RDX  molecules,  although  insight  into  the 
nature  of  the  interactions  was  obtained  through  calculation  of  the  DMNA-DMNA  interactions. 
Thus,  a  functional  form  for  the  RDX-RDX  interactions  was  assumed,  and  the  resulting  model 
subjected  to  a  variety  of  molecular  simulations  in  order  to  assess  its  quality.  The  intermolecular 
interactions  between  the  molecules  of  the  crystal  are  approximated  using  superpositions  of 
pairwise  Buckingham  (6-exp)  (repulsion  and  dispersion)  and  Coulombic  (C)  potentials  of  the 
form: 

Vapfr>  =  Aa|?exP(-Bafir)-cal!/r6"  « 

and 


Vc  /_>,  q«qP 
V“P(r)-4^0r 


(10) 


where  r  is  the  interatomic  distance  between  atoms  a  and  p,  qa  and  qp  are  the  electrostatic  charges 
on  the  atoms,  and  eo  is  the  dielectric  permittivity  constant  of  space. 


The  parameters  for  the  6-exp  potential  in  Eq.  (9)  were  adjusted  to  reproduce  the 
experimental  crystal  structure  of  RDX  at  ambient  conditions,  and  the  partial  charges  used  in  the 
Coulombic  interactions  [Eq.  (10)]  were  those  that  reproduced  the  quantum  mechanically  derived 
electrostatic  potential  surrounding  the  isolated  molecule.  The  resulting  model  proved  to  be 
extremely  accurate  in  its  description  of  RDX  at  a  large  range  of  temperatures  and  pressures  [20], 
and  also  proved  to  be  transferable;  that  is,  it  was  able  to  describe  a  large  number  of  nitramine 
[21-23]  and  non-nitramine  energetic  molecular  crystals[24-25].  The  results  will  be  described 
hereafter. 


The  first  model  used  to  describe  RDX-CO2  interactions  assumes  a  potential  function  V(R) 
that  has  a  Lennard-Jones  form,  where  R  is  the  separation  between  molecular  centers-of-mass.  Its 
parameters  were  selected  according  to  the  following  procedure.  A  rotationally-averaged 
interaction  potential  V’(R)  is  calculated: 

V(R)-<E£vjk(Rjl£)>.  (11) 

j  k 

Vjk  is  the  interaction  potential  between  j-th  atom  of  RDX  and  the  k-th  atom  of  the  CO2  molecule, 
where  the  atom  pairs  are  separated  by  a  distance  Rjk  and  the  molecules  are  arranged  relative  to 
one  another  in  an  arbitrary  orientation.  The  potential  term  V’(R)  represents  the  average  over  a 
large  number  (~  5xl05)  of  random  relative  orientations  of  RDX  and  the  CO2  molecule  with  the 
centers  of  mass  separated  by  the  distance  Rj.  The  atom-atom  interactions  in  these  calculations 
are  represented  as  Buckingham  potentials  of  the  form  described  in  Eq.  (9).  The  parameters  Ajk, 
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Bjk,  and  Cjk  in  the  above  equation  are  those  used  to  predict  the  crystal  packing  of  RDX  [20].  For 
the  mixed  atomic  interactions  we  have  used  the  arithmetic  mean  rule  for  B  and  the  geometric 
mean  rule  for  the  A  and  C  parameters. 

We  originally  selected  Lennard-Jones  parameters  that  reproduced  the  location  of  the 
minimum  and  well  depth  of  the  result  of  Eq.  (11);  however,  solubility  predictions  using  this 
model  were  too  low  by  several  orders  of  magnitude  compared  to  the  experimental  extraction 
results.  We  found  empirically  that  a  Lennard-Jones  potential  whose  minimum  and  well  depth 
were  shifted  from  those  of  Vj’(Rj)  by  factors  of  0.97  and  1.20,  respectively,  significantly 
improved  the  agreement  between  the  calculated  and  experimental  solubility  values.  The  RDX- 
CO2  interactions  are  described  by  a  L-J  potential  with  parameters  <7=5.282  A  and  £=0.792 
kcal/mol .  The  CO2-CO2  interactions  are  described  by  the  Moller-Fischer  potential  [26]. 

The  following  assumptions  were  made  for  the  description  of  interactions  involving  the 
modifier  molecules.  In  this  model  a  modifier  molecule  is  considered  to  consist  of  an  ensemble  of 
m  groups.  For  example,  for  DMSO  there  are  m=3  such  groups,  i.e.  CH3,  CH3  and  SO.  Similarly, 
for  CH3OH  and  CH3CN  there  are  two  groups  (m=2),  respectively  (CH3,  OH)  and  (CH3 ,  CN). 

The  interaction  potential  between  RDX  and  a  modifier  molecule  is  taken  as  a  sum  of  the 
following  terms: 


Vrdx-m_Vdd+  ^  Vj(Rj),  (12) 

where  Vdd  is  a  dipole-dipole  interaction  term,  Vi(Ri)  gives  the  interaction  between  the  RDX 
molecule  and  the  i-th  group  of  the  modifier  molecule  and  Rj  is  the  distance  between  the  centers 
of  mass  of  RDX  and  the  group.  The  summation  is  taken  over  all  m  groups  of  the  modifier 
molecule. 

The  potential  Vi(Ri)  has  a  Lennard-Jones  form,  and  its  parameters  were  selected 
according  to  the  procedure  developed  to  describe  the  RDX-CO2  interaction,  described  above.  In 
this  procedure,  a  rotationally-averaged  interaction  potential  Vj’(Ri)  is  calculated: 

Vl'<Rl)  =  <LEVjk(Rjk)>„  (13) 

j  k 

Vjk  is  the  interaction  potential  between  j-th  atom  of  RDX  and  the  k-th  atom  of  the  group  placed 
at  the  distance  Rjk  with  an  arbitrary  orientation.  The  potential  term  Vj’(Ri)  represents  the  average 
over  a  large  number  (~  5x10s)  of  random  relative  orientations  of  RDX  and  the  group  with  the 
centers  of  mass  separated  by  the  distance  Rj.  The  atom-atom  interactions  in  these  calculations 
are  represented  as  Buckingham  potentials  of  the  form  in  Eq.  (9),  and  the  parameters  for  these  are 
are  those  used  to  predict  the  crystal  packing  of  RDX  [20].  For  the  mixed  atomic  interactions  we 
have  used  the  arithmetic  mean  rule  for  B  and  the  geometric  mean  rule  for  the  A  and  C 
parameters. 

For  consistency,  we  have  used  the  same  scaling  factors  to  obtain  Lennard-Jones 
parameters  to  describe  the  RDX  interactions  with  the  groups  on  the  modifier  molecules  that  were 
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used  for  the  RDX  interactions  with  CO2.  The  same  parameters  were  used  to  describe  the  CH3 
groups  on  DMSO,  CH3CN  and  CH3OH;  however,  these  were  obtained  using  the  results  of  Eq. 
(13)  for  the  DMSO  interactions  with  RDX. 

For  the  dipole-dipole  interaction  term,  Vdd,  we  consider  that  both  RDX  and  the  modifier 
can  be  considered  as  point  dipoles.  In  this  case,  the  dipole-dipole  interaction  is  given  by  an 
angle-dependent  term  [9]  proportional  to  (pa.pt/R3) ,  where  pa  and  p.b  denote  the  dipole 
moments  of  RDX  and  the  modifier,  respectively,  and  R  gives  the  distance  between  the  centers 
of  mass  of  RDX  and  the  modifier. 

In  describing  the  C02-modifier  and  modifier-modifier  interactions,  the  molecules  are 
considered  to  be  point  masses.  The  interactions  are  described  by  the  Lennard-Jones  function, 
with  the  parameters  obtained  in  the  same  manner  as  those  obtained  for  RDX-CO2,  including 
using  the  same  scaling  factors  (0.97  for  cr  and  1.20  for  e).  Thus  this  model  considers  a  non¬ 
isotropic  interaction  for  each  RDX-modifier  pair  and  isotropic  interactions  for  C02-modifer  and 
modifier-modifier  pairs. 

A  second  model  interaction  potential  was  developed  to  describe  the  RDX-CO2 
interactions  to  investigate  the  effect  of  assuming  a  multiple-site  model  of  this  interaction  (as  used 
for  the  modifiers).  In  this  model,  one  center  is  located  on  the  CO2  molecule  and  the  RDX 
molecule  has  seven  sites,  each  located  at  the  center  of  mass  of  its  group.  The  groups  consist  of 
the  ring,  the  three  NO2  groups,  and  the  three  hydrogen-atom  pairs.  Parameters  were  determined 
in  the  manner  described  for  the  RDX-modifier  interactions;  however,  the  scaling  factors  differ. 

In  this  case,  the  scaling  factors  for  e  and  a  are  1 .55  and  0.88,  respectively.  This  second  model 
will  be  denoted  as  the  “multisite”  model  (compared  to  the  previously-described  monosite  model). 
Simulations  using  this  model  of  interactions  of  RDX-CO2  were  not  used  in  simulations  in  which 
polar  modifiers  were  introduced  into  the  solvent. 

TNT/CO2  systems 

Two  models  to  describe  the  interactions  between  TNT  and  CO2  were  developed  in  a 
manner  similar  to  that  for  RDX  and  CO2:  The  first  assuming  a  single  interaction  center  on  the 
TNT  and  CO2  molecules,  respectively,  and  the  second  assumes  seven  sites  on  the  TNT  molecule 
(the  ring,  the  methyl  group,  the  three  NO2  groups,  and  the  two  hydrogen  atoms)  and  one  on  the 
CO2  molecule.  The  monosite  interaction  potential  for  TNT  with  CO2  requires  that  the  scaling 
parameters  be  different  from  those  used  in  the  RDX-CO2  interaction  in  order  to  obtain  good 
agreement  with  the  experimental  solubility  data. 

Results 


Numerous  computer  experiments  were  performed  to  assess  the  quality  of  the  interaction 
potentials  that  were  developed  and  to  establish  the  dependencies  on  the  enhancement  of 
solubility  due  to  the  composition  and  nature  of  the  solvent. 

Crystalline  RDX 
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This  model  was  developed  to  achieve  modest  and  narrow  goals,  i.e.  to  study  the  nitramine 
explosive  RDX  (l,3,5-hexahydro-l,3,5,-x-triazine)  [20].  Isothermal-isobaric  molecular  dynamics 
simulations  (NPT-MD)  and  molecular  packing  calculations  (MP)  were  used  to  assess  the 
interaction  potential  through  the  prediction  of  the  crystal  structure  of  RDX  at  ambient  conditions. 
The  results  showed  that  this  model  reproduced  the  crystallographic  parameters  of  RDX  to  within 
2%  of  experiment  in  the  cell  dimensions,  with  little  rotational  or  translational  disorder  within  the 
unit  cell  [20].  Also,  we  found  that  this  interaction  potential  could  also  describe  different 
polymorphic  phases  of  two  other  nitramine  crystals:  the  polycyclic  nitramine  2,4,6,8,10,12- 
Hexanitrohexaazaisowurtzitane  (HNIW,  or  CL-20)  [21]  and  the  monocyclic  nitramine 
Octahydro-l,3,5,7-Tetranitro-l,3,5,7-tetraazacyclooctane  (HMX)  [22],  Molecular  simulations 
for  these  crystals  predicted  cell  parameters  within  a  few  percent  of  experiment,  and  little 
translational  or  rotational  disorder  of  the  molecules.  These  successes  called  for  further 
investigation  to  determine  the  limits  of  the  transferability  of  this  interaction  potential  to  other 
energetic  crystals.  We  performed  molecular  packing  (MP)  calculations  of  30  nitramine  crystals 
[23],  For  most  of  these  crystals,  the  predicted  structural  lattice  parameters  deviate  by  less  than 
5%  from  experiment.  Also,  for  most  of  the  crystals  there  are  small  rotations  and  practically  no 
translations  for  the  molecules  in  the  asymmetric  unit  cell.  The  interaction  potential  was  further 
assessed  in  molecular  packing  calculations  of  5 1  crystals  containing  non-nitramine  molecules 
with  functional  groups  common  to  energetic  materials  [24].  MP  calculations  using  this 
interaction  potential  reproduced  the  crystal  structures  to  within  5%  of  experiment  for  these  5 1 
non-nitramine  systems,  including  the  explosives  PETN,  nitromethane,  TATB,  TNT,  TNAZ  and 
several  nitrocubane  derivatives.  Also,  MD  simulations  of  RDX,  HMX,  and  HNIW  under 
hydrostatic  compression  produced  results  that  were  in  good  agreement  with  experiment  over  the 
range  of  pressures  investigated  experimentally  [25].  MD  simulations  of  PETN  under  hydrostatic 
compression  were  in  acceptable  agreement  with  experiment  up  to  5  GPa;  beyond  that,  the 
disagreement  in  predictions  with  experiment  was  attributed  to  the  inadequacy  of  the  rigid-body 
approximations  imposed  in  these  simulations  when  applied  to  floppy  molecules.  The 
significance  of  this  effort  lies  in  the  ability  of  the  model  to  predict  properties  for  a  wide  range  of 
energetic  crystals,  should  such  be  desired  by  the  user. 

RDX  solubility  in  pure  SF  CO2 

Calculated  values  of  the  solubility  of  RDX  in  pure  SF  CO2  using  the  monosite  and 
multisite  interaction  models  for  the  RDX-CO2  interactions  are  given  in  Table  1 .  Experimental 
data  of  Morris  [1]  are  also  given  for  comparison.  The  general  trend  in  the  experimental  and 
computed  data  at  a  given  temperature  is  that  the  solubility  increases  with  pressure.  At  P  >  13.8 
MPa,  the  solubility  exhibits  a  maximum  as  a  function  of  temperature.  The  results  for  P  >  13.8 
MPa  suggest  that  there  may  be  a  maximum  in  the  solubility  as  a  function  of  temperature  beyond 
the  range  of  temperatures  considered  in  this  study.  Qualitatively,  this  behavior  is  explained  by 
Eq.  (6).  The  solubility  is  both  temperature-  and  density-dependent  directly  and  indirectly 
[through  its  dependence  on  p*  (which  is  temperature  dependent)  and  pr  (which  is  density 
dependent)] .  As  the  temperature  increases,  there  is  an  increase  in  p*,  which,  according  to  Eq. 

(6),  increases  the  value  of  the  solubility.  At  the  same  time,  the  density  decreases,  resulting  in  a 
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decrease  in  the  value  of  -p.r  and  the  solubility.  Thus,  there  are  competing  effects  with  increasing 
temperature  that  could  result  in  a  maximum  on  the  solubility  vs.  temperature  curve. 


A  comparison  of  the  experimental  and  computed  solubilities  in  Table  1  shows  that  within 
the  limits  of  the  uncertainties  in  the  computed  values,  there  is  very  good  agreement.  Such 
agreement  in  the  solubility  data  over  the  large  range  of  0.001  to  0.254  mg  RDX/g  CO2  is  notable, 
since  a  small  error  in  the  residual  chemical  potential  can  lead  to  a  large  error  in  the  solubility. 

The  computed  results  show  the  same  trends  in  the  solubility  as  functions  of  temperature  and 
pressure  as  observed  by  experiment.  For  the  majority  of  the  points,  the  agreement  between 
experiment  and  theory  is  within  a  factor  of  1.5.  The  worst  agreement  corresponds  to  the 
solubility  at  points  (338  K,  10.4  MPa)  and  (323  K,  10.4  MPa)  where  the  experimental  solubility 
values  are  very  low.  Computed  values  by  both  models  are  lower  than  the  experimental  values. 

Table  1.  Solubility  of  RDX  in  CO2  as  a  function  of  pressure  and  temperature. 


Temp.  (K) 

Solubility  (mg/g) 

Multisite  Monosite  Expt.1,2 

353.0 

48.3 

0.873 

0.270 

0.257 

0.254 

41.4 

0.840 

0.197 

0.237 

27.6 

0.732 

0.081 

0.114 

13.8 

0.373 

0.004 

338.0 

48.3 

0.919 

0.173 

41.4 

0.889 

0.113 

0.173 

27.6 

0.806 

0.068 

| 

0.076 

13.8 

0.006 

- 

0.009 

10.4 

0.278 

0.0001 

0.0001 

■sum 

323.0 

48.3 

0.963 

0.116 

0.193 

— 

41.4 

0.938 

■MB 

0.092 

■E2B 

27.6 

0.871 

■E5B 

- 

m-zm 

13.8 

0.720 

0.011 

0.021 

0.013 

10.4 

0.443 

0.0005 

0.0003 

0.003 

308.0 

48.3 

HHESSIB 

0.005 

fffESB 

41.4 

■EBB 

0.026 

27.6 

0.931 

BKIilHl 

0.040 

■tfilcgM 

13.8 

0.823 

■BUB 

- 

0.013 

0.784 

■KB 

0.015 

0.008 

303.0 

48.3 

1.022 

■i  1 

0.007 

■EeEB 

41.4 

BBEZ2SB 

0.025 

0.006 

■EiB 

27.6 

BKEZSBB 

0.023 

- 

0.032 

13.8 

0.867 

0.011 

- 

0.013 

10.4 

0.824 

0.005 

0.009 

0.007 

A  comparison  of  the  results  of  the  two  models  shows  that  the  results  of  the  monosite 
model  deviate  from  experiment  and  results  from  the  multisite  model  at  low  temperatures  and 
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high  pressures  corresponding  to  high  densities  of  the  solvent.  The  disagreement  at  high  densities 
is  attributed  to  differences  in  the  interaction  potentials  for  intermolecular  separations  smaller  than 
6  A.  For  the  two  models,  the  interaction  potentials  between  molecular  centers  are  identical  for 
separations  greater  than  these;  for  smaller  separations,  the  differences  are  significant. 

RDX  solubility  in  polar-modified  SF  CO2 

We  have  attempted  to  identify  the  molecular  properties  of  representative  co-solvents  that 
affect  the  solubility  enhancement  of  RDX  in  polar  modified  SF  CO2.  The  modifiers,  DMSO, 
CH3CN  and  CH3OH,  were  chosen  to  represent  very  effective,  moderately  effective,  and  non- 
effective  co-solvents  for  solubility  enhancement  of  RDX  in  SF  CO2.  The  most  obvious 
distinctions  between  these  modifiers  are  the  molecular  sizes  and  shapes  and  the  magnitude  of  the 
dipole  moments.  The  dipole  moments  of  DMSO  and  CH3CN  (3.96  D  and  3.92  D,  respectively) 
are  very  similar,  but  nearly  twice  as  large  as  that  of  CH3OH  (1.70  D)  [27]. 

(a)  Modifier  dependence 

Table  2  gives  the  computed  solubility  results  for  all  three  modifiers  at  T  =  323  K  and 
P=41.4  MPa.  The  solubility  enhancement  factor  (SEF)  is  the  ratio  of  the  solubility  of  RDX  in  the 
polar-modified  CO2  to  the  value  for  RDX  in  pure  CO2  (0.0916  mg  RDX/g  CO2  [10]).  The  SEFs 
reported  in  the  table  show  the  proper  ranking  in  the  effectiveness  of  the  polar  modifiers,  i.e. 
DMSO  >  CH3CN  >  CH3OH  compared  to  the  experimental  EEFs.  The  results  also  qualitatively 
predict  the  degree  of  enhancement  as  a  function  of  modifier.  It  is  expected  that  the  experimental 
EEFs  would  be  smaller  than  the  true  SEFs  due  to  the  experimental  procedure  and 
multicomponent  nature  of  the  sample  used  in  the  experiments.  Thus  exact  agreement  between 
our  calculations  and  experiment  would  not  be  expected. 

(b)  Concentration  dependence 

Table  2  also  shows  the  variation  in  the  solubility  enhancement  factor  as  a  function  of 
mole  concentration  for  DMSO  and  CH3CN.  Figure  1  illustrates  this  information  for  DMSO-CO2 
solvents  and  provides  a  comparison  of  the  degree  of  enhancement  predicted  by  theory  (upper 
frame)  with  that  observed  in  the  experimental  extractions  (lower  frame).  Although  the 
magnitudes  of  the  enhancements  differ  between  the  predictions  and  experiment,  the  behavior  of 
the  enhancement  with  increasing  modifier  concentration  is  the  same.  The  enhancement  factor 
increases  from  14.3  at  2%  DMSO  mole  concentration  to  700.6  at  6%  DMSO  mole  concentration. 
Figure  (1)  shows  that  with  the  increase  in  concentration  there  is  a  rapid  (exponential)  increase  in 
the  enhancement  factor  followed  by  a  slower  increase,  suggesting  that  the  limit  of  saturation  is 
being  approached. 

We  also  performed  solubility  calculations  in  which  the  dipole-dipole  interaction  was  set 
to  zero  for  all  concentrations  of  DMSO  reported  in  Table  2.  The  calculated  enhancement  factors 
were  all  close  to  unity,  and  ranged  between  0.65  and  1.57.  This  indicates  that  the  increase  in 
enhancement  factor  with  increasing  concentration  is  mainly  due  to  the  dipole-dipole  interaction 
term.  For  higher  concentrations  of  the  modifier,  the  model  requires  a  further  modification  to 
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account  for  the  rotation  of  the  modifier  molecules  and  a  correlation  in  the  directions  of  the  dipole 
moments  of  the  neighboring  molecules. 

Table  2.  Solubility  and  solubility  enhancement  factor  (SEF)  of  RDX  due  to  polar  modifiers 


Mole 

Concentration 

Partial 
Density  of 
C02  (g/cc) 

Theoretical 

Experimental 

[1,21 

SEF 

Amount 
extracted 
(mg  RDX/g 
CO2) 

EEF 

DMSO 

0.0 

0.937 

1.0 

0.01  l,a 
0.08b 

0.899 

14.3 

3.38 

0.876 

85.9 

0.07 6a 

69a 

3.5 

mamm 

41b 

3.99 

0.867 

176.2 

5.00 

6.9b 

~86* 

6.00 

0.837 

6.30 

9.5b 

119* 

CH3CN 

0.0 

0.937 

1.0 

HHE33HI 

3.8 

12b 

3.99 

0.835 

24.4 

25a 

6.00 

0.793 

174.7 

6.1 

1.3b 

16b 

CH3OH 

0.937 

1.0 

5.54 

0.834 

2.7 

1.4a 

a.  Static-dynamic  extraction  on  bulk  M43  Propellant  (Ref.  1) 

b.  Dynamic  extraction  on  ground  M43  Propellant  (Ref.  2) 

(c)  Dipole  moment  dependence 

We  have  investigated  the  effect  of  changing  pb  while  keeping  all  other  parameters 
(including  the  L-J  interaction  parameters)  unchanged  in  order  to  establish  the  effect  of  dipole- 
dipole  interactions  between  RDX  and  a  polar  modifier.  We  performed  calculations  corresponding 
to  experimental  conditions  (T=323  K,  P=41.4  MPa)  for  DMSO  (3.38%),  CH3CN  (3.99%)  and 
CH3OH  (5.54%),  but  the  dipole  moment  was  set  to  zero  for  each  of  the  modifiers.  The  other 
parameters  of  the  potential  were  not  changed  in  the  simulation.  This  has  the  effect  of  including 
molecular  size  and  shape,  but  it  precludes  electrostatic  interactions.  The  DMSO,  CH3CN  and 
CH3OH  systems  have  enhancement  factors  of  1 .4,  1 .0  and  0.7,  respectively,  under  these 
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conditions.  This  suggests  that  the  presence  of  a  dipole  moment  has  a  significant  effect  on 
solubility  enhancement. 

We  next  investigated  the  variation  in  the  solubility  enhancement  factor  as  a  function  of  pb 
for  CO2-DMSO  (3.38%)  at  T=323  K  and  P=4 1 .4  MPa.  The  results  are  given  in  Table  3  and 
shown  in  Fig.  2.  These  data  suggest  that  the  enhancement  factor  is  very  sensitive  to  the 
magnitude  of  the  dipole  moment  of  the  modifier.  Figure  2  shows  that  the  increase  in  the 
solubility  with  the  increase  in  the  dipole  moment  of  the  modifier  is  nearly  exponential.  The 
enhancement  factor  for  a  null  dipole  moment  is  1 .4  and  is  870.3  for  a  dipole  moment  equal  to  5.0 
D.  This  suggests  that  the  low  value  of  the  enhancement  factor  for  the  CH3OH/CO2  system  is  due 
mainly  to  the  small  dipole  moment  (1 .70  D). 

Table  3.  Solubility  Enhancement  Factor  (SEF)  as  function  of  modifier  dipole  moment 
(3.38%  concentration  of  modifier  in  CO2,  T=323  K,  P=41.4  MPa).  All  other  modifier 
parameters  correspond  to  DMSO. 


pb  (Debye) 

SEF 

0.0 

1.4 

1.0 

1.8 

2.0 

4.0 

3.0 

15.1 

3.96 

85.9 

5.0 

870.3 

(d)  Comparison  of  the  Enhancement  Factors  of  DMSO  and  CH3CN 

The  preceding  discussions  indicate  that  the  magnitude  of  the  dipole  moment  is  a 
significant  factor  in  solubility  enhancement,  however,  it  does  not  account  for  the  total  effect.  This 
is  evident  upon  comparison  of  the  performance  of  DMSO  and  CH3CN  as  polar  modifiers. 

DMSO  and  CH3CN  have  nearly  the  same  values  of  dipole  moment.  If  the  solubility 
enhancement  factor  is  a  function  of  only  the  dipole  moment  of  the  modifers,  then  CH3CN  should 
have  a  solubility  enhancement  factor  that  is  only  5%  smaller  than  DMSO,  according  to  the  results 
shown  in  Table  3  and  Fig.  2.  However,  the  enhancement  factors  at  3.99%  modifier  concentration 
for  CH3CN  (24.5)  is  86%  smaller  than  that  of  DMSO  (176.2).  Clearly,  there  are  additional 
factors  other  than  dipole  moment  which  affect  the  solubility.  One  factor  could  be  the  effect  on 
the  solvent  structure  due  to  the  modifier.  We  will  investigate  the  solvent  structure  by  analyzing 
energetically  favored  configurations  of  the  test  particle-solvent  system  and  by  calculating  pair 
distribution  functions  for  the  solvent  molecules. 

Within  these  simulations  we  know  that  the  largest  contributors  to  the  value  of  pr  are  those 
configurations  for  which  the  potential  'P  experienced  by  the  test  particle  is  attractive  [Eq.  (7)]. 
Likewise,  we  know  the  configurations  that  do  not  have  a  significant  contribution  to  the  value  of 
qr  are  those  for  which  TX),  which  are  typically  due  to  repulsive  interactions  by  close  approach  of 
solvent  molecules  to  the  test  particle.  We  have  extracted  from  the  total  those  configurations  of 
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the  RDX/C02-modifier  solution  that  have  values  of  TO.  The  number  of  such  configurations  is 
small;  for  example,  in  a  typical  run  we  found  that  the  number  of  configurations  for  which  T  is 
attractive  (T  <  0 )  is  34  million  out  of  a  total  250  billion  insertions. 

We  next  attempted  to  quantify  the  number  of  modifier  molecules  in  the  neighborhood  of 
RDX  for  the  energetically-favored  configurations.  For  this  purpose,  we  considered  a  modifier 
molecule  to  be  a  neighbor  if  it  is  within  7.0  A  distance  of  the  RDX  molecule.  Figure  3  shows 
the  average  number  of  DMSO  neighbors  of  a  RDX  molecule  as  a  function  of  T  for  fluid 
configurations  at  T  =  323  K,  P=41.4  MPa,  and  DMSO  concentration  6%.  The  figure  shows  that 
lowest-energy  configurations  are  those  for  which  there  are  more  than  one  modifier  in  the 
neighborhood  of  the  RDX  particle.  The  behavior  of  the  curve  suggests  a  preference  of  RDX  for 
occupying  a  site  that  has  three  or  more  modifier  molecules  as  its  neighbors. 

The  pair  distribution  function  g(R)  for  CO2-CO2,  DMSO-DMSO  and  DMSO-CO2  for 
CO2  with  3.38%  DMSO  at  P=41.4  MPa  and  T=323  K  were  calculated,  and  shown  in  the  upper 
frame  of  Fig.  4.  The  corresponding  curves  for  CO2-CH3CN  (3.99%)  and  CO2-CH3OH  (5.54%) 
at  the  same  P  and  T  are  shown  in  the  middle  and  lower  frames  of  Fig.  4,  respectively.  The 
leading  peak  for  CO2-CO2  in  all  the  three  curves  is  located  at  4.05+/-  0.05  A.  The  location  of  the 
leading  peaks  for  the  C02-modifier  and  modifier-modifier  pair  distribution  functions  varies  with 
the  modifier.  The  location  of  the  leading  peak  is  5  to  9%  larger  than  the  corresponding  value  of 
cr  parameter  of  L-J  interaction.  The  locations  of  these  peaks  are  insensitive  to  the  modifier 
concentration,  temperature  and  pressure  in  the  range  of  these  studies.  These  figures  show  that 
CH3CN-CH3CN  intermolecular  separations  are  about  12%  smaller  than  the  DMSO-DMSO 
separations.  This  is  consistent  with  the  10%  higher  value  of  the  sigma  parameter  for  the  DMSO- 
DMSO  interactions  than  that  for  CH3CN-CH3CN.  According  to  the  arithmetic  mean  rule  for  the 
sigma  parameter,  a  10%  difference  in  the  modifier-modifier  interactions  would  result  in  a  5% 
difference  in  the  effective  sizes  of  RDX-CH3CN  and  RDX-DMSO. 

The  difference  in  the  RDX-CH3CN  and  RDX-DMSO  interactions  can  also  be  understood 
from  Fig.  5.  Figure  5  shows  the  interaction  potential  given  by  Eq.  (12)  between  RDX  and  a 
modifier  molecule  as  a  function  of  the  distance  R  between  their  centers  of  mass  in  the  most 
favorable  orientation.  The  most  favorable  orientation  is  defined  here  as  the  orientation  of  the 
modifier  molecule  relative  to  RDX  such  that  the  directions  of  the  dipole  moment  vectors  are 
antiparallel  and  the  potential  given  by  Eq.  (12)  is  the  minimum  of  all  values  obtained  for  1 05 
random  relative  orientations  of  the  pair  for  a  given  value  of  R. 

Figure  5  shows  that  the  value  of  the  well-depth  of  such  a  potential,  denoted  as  V(min.),  is 
6.655  kcal/mole  at  R  =5.38  A  for  the  DMSO-RDX  system.  The  well-depth  of  V(min.)  for  the 
CH3CN-RDX  system  is  6.385  kcal/mole  at  R=5.14  A.  As  expected,  the  effective  size  parameter 
for  the  RDX-DMSO  interaction  is  larger  by  about  5%.  Also,  the  effective  interaction  strength  is 
larger  for  the  RDX-DMSO  pair  by  ~  0.27  kcal/mole.  This  difference  in  the  interaction  strength  is 
one  of  the  important  factors  for  the  higher  solubility  of  RDX  in  DMSO.  Additionally,  the  size 
difference  could  explain  why  RDX  is  less  soluble  in  CH3CN.  We  illustrate  this  proposal  using 
Fig.  6.  The  three  circles  in  the  left-hand  portion  of  this  figure  depict  three  DMSO  molecules  that 
surround  the  RDX  molecule.  This  is  a  simple  two-dimensional  representation  of  a  low-energy 
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fluid  configuration  determined  earlier  (energetically  favored  configurations  are  those  in  which 
three  or  more  modifiers  are  in  the  neighborhood  of  RDX).  The  four  molecules  may  or  may  not 
lie  in  the  same  plane.  The  solid  lines  connecting  the  modifier  molecules  represent  the  most 
probable  distances  between  DMSO-DMSO  molecules. 

The  right-hand-side  of  the  figure  represents  the  CH3CN-RDX  analog.  The  three  circles  at 
each  apex  of  the  triangle  connected  by  solid  lines  represent  separations  that  are  12%  smaller  than 
those  of  the  DMSO  system.  This  is  the  same  difference  in  the  molecular  separations  that  were 
determined  from  the  radial  distribution  functions  in  Fig.  5.  As  discussed  in  preceding 
paragraphs,  however,  we  estimated  that  the  energetically-favored  CH3CN-RDX  distances  are 
only  5%  smaller  than  the  DMSO-RDX  distances,  not  12%  smaller.  Therefore  an  RDX  molecule 
inserted  into  a  region  in  which  it  was  surrounded  by  three  CH3CN  molecules  would  have  a 
higher  probability  of  experiencing  a  repulsive  interaction  than  if  it  was  surrounded  by  three 
DMSO  molecules  in  a  CO2-DMSO  solvent.  The  increased  number  of  repulsive  states  sampled  in 
a  simulation  results  in  a  lower  solubility.  Repulsive  interactions  can  be  avoided  if  the  test 
particle  is  inserted  at  a  location  that  is  a  sufficient  distance  from  the  three  modifier  molecules. 
Such  an  arrangement  may  provide  a  non-negligible  contribution  to  the  average  in  Eq.  (7). 
However,  the  RDX  test  particle  may  experience  repulsions  from  other  CH3CN  molecules  in  the 
solution.  Therefore,  the  test  particle  samples  more  energetically  unfavorable  sites  in  a 
CO2/CH3CN  solution  than  in  a  CO2/DMSO  solution.  This  explanation  suggests  that  in  an 
experiment,  the  number  of  arrangements  of  RDX  surrounded  by  a  modifiers  at  molecular 
separations  that  give  the  optimal  attraction  for  both  modifier-modifier  and  RDX-modifier 
interactions  is  smaller  for  CH3CN  than  DMSO,  hence  the  lower  solubility.  To  summarize: 
Although  the  dipole  moments  of  CH3CN  and  DMSO  are  nearly  same,  their  contributions  to 
solubility  enhancements  of  RDX  at  the  P,  T,  and  concentration  of  our  studies  are  different 
because  of  the  difference  in  (1)  RDX-modifier  interaction  strength  and  (2)  the  size  of  the 
modifiers. 

(e)  Dependence  on  temperature  and  pressure 

The  computed  solubility  of  RDX  in  CO2-DMSO  (3.38%)  as  a  function  of  temperature 
and  pressure  for  a  few  typical  values  of  P  and  T  are  reported  in  Table  4.  The  solubility 
enhancement  factors  relative  to  the  computed  solubility  values  in  pure  CO2  [10]  have  also  been 
reported  in  the  table.  The  results  show  that  the  solubility  increases  with  the  pressure.  A 
comparison  of  the  solubility  of  RDX  in  CO2  with  and  without  [10]  modifier  shows  that  the 
solubility  of  RDX  in  polar-modified  CO2  is  less  sensitive  to  temperature. 
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Table  4.  Solubility  and  Solubility  Enhancement  Factor  (SEF)  of  RDX  in  CO2  with  3.38% 
mole  concentration  of  DMSO  as  a  function  of  temperature  and  pressure 


T  (K) 

P  (MPa) 

BEE358 

IMiMH 

SEF 

323.0 

13.8 

0.698 

0.833 

39.4 

41.4 

0.876 

7.868 

85.9 

48.3 

0.900 

55.2 

338.0 

41.4 

0.836 

9.046 

47.1 

48.3 

0.863 

12.315 

74.2 

353.0 

41.4 

0.794 

7.419 

41.7 

48.3 

0.823 

9.462 

37.0 

(f)  Other  models 

We  also  calculated  the  solubility  at  the  temperature  and  pressure  given  in  Table  2  for  all 
three  modifiers  using  models  that  treat  each  modifier  molecule  as  a  point  mass,  i.e.,  taking  m=l 
in  Eq.  (12).  In  this  model,  all  potential  energy  features,  including  the  random  selection  of 
directions  of  dipole  moment  vectors,  are  determined  by  the  same  procedure  described  earlier. 

For  purposes  of  discussion,  this  model  will  be  denoted  as  “Model  B”,  whereas  the  model 
described  in  the  earlier  sections  will  be  denoted  as  “Model  A”.  The  calculated  enhancement 
factors  for  Model  B  are  1.7, 6.2  and  0.7,  for  DMSO(3.38%),  CH3CN(3.99%),  and 
CH30H(5.54%).  The  inability  of  this  model  to  qualitatively  reproduce  the  behavior  of 
experimental  EEF  may  be  due  to  either  the  inappropriate  selection  of  potential  scaling  parameters 
(0.97,  1.20  for  a  and  e,  respectively)  or  to  the  failure  of  the  description  of  RDX-modifier 
interaction  by  the  sum  of  a  single  L-J  potential  term  and  a  dipole-dipole  interaction  term. 

A  third  model,  denoted  “Model  C”,  assumes  modifier-modifier  and  C02-modifier 
interactions  that  are  different  from  those  of  Model  A.  The  CCVmodifier  and  modifier-modifier 
interactions  in  Model  C  require  that  each  modifier  have  m  interaction  sites,  similar  to  the 
requirements  for  the  RDX-modifier  interactions  as  described  earlier.  Use  of  this  model  in  the 
NPT-MC  generation  of  solvent  configurations  does  not  allow  random  assignment  of  the  dipole 
moment  of  the  modifier  molecule,  as  described  in  simulations  using  Model  A. 

The  enhancement  factors  calculated  using  Model  C  are  found  to  be  very  low.  The 
predicted  enhancement  factor  for  CO2  with  3.38%  DMSO,  T=323  K,  P=41.4  MPa  is  0.7.  We 
have  repeated  calculations  using  Model  C  for  different  sets  of  potential  scaling  parameters 
(ranging  from  0.97  to  1 .35  for  cr  and  from  0.80  to  1 .20  for  s)  for  the  C02-modifier  and  modifier- 
modifier  potentials  such  that  the  densities  of  the  fluid  using  these  different  potential  scaling 
parameters  are  almost  the  same  as  that  given  in  Table  2  by  Model  A.  However,  the  enhancement 
factors  were  very  low  in  all  cases.  Further,  solubilities  predicted  using  Model  C  show  a  weak 
dependence  on  dipole  moment,  e.g.,  the  enhancement  factor  varies  from  0.4  to  1 .9  as  the  dipole 
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moment  of  the  modifier  is  varied  from  0.0  to  6.0  D,  while  keeping  all  other  parameters  the  same 
as  that  for  DMSO. 


We  think  the  differences  in  behavior  between  Model  C  and  Model  A  can  be  attributed  to 
the  following:  In  Model  C,  the  configurations  corresponding  to  the  lowest  enthalpy  values  are 
those  in  which  the  dipole  moments  of  adjacent  modifier  molecules  are  aligned  in  an  antiparallel 
orientation.  In  such  a  configuration,  the  dipole-dipole  interaction  between  the  test  particle 
(RDX)  and  a  modifier  neighbor  would  tend  to  cancel  the  dipole-dipole  interaction  between  the 
RDX  and  the  adjacent  modifier  neighbor  whose  dipole  is  in  the  antiparallel  orientation.  Such  a 
cancellation  would  lead  to  the  lower  contribution  to  the  solubility.  In  Model  A,  on  the  other 
hand,  the  orientations  of  the  dipole  moments  are  chosen  randomly.  Hence  the  random  selection 
of  the  direction  of  the  dipole  moments  of  two  modifier  neighbors  may,  in  many  cases,  have 
favorable  orientations  for  which  the  RDX-modifier  dipole-dipole  interactions  are  not  cancelled. 

At  a  first  glance,  Model  C  appears  to  be  more  physically  representative  of  the  actual  fluid 
than  Model  A.  However,  if  we  consider  continuous  thermal  rotations  of  the  molecules  in  the  real 
fluid,  then  we  may  find  that  the  modifier  molecules  rotate  with  different  frequencies  and  as  such 
the  relative  orientations  of  dipole  moment  vectors  of  the  two  modifier  molecules  may  frequently 
become  parallel  and  antiparallel.  Since  the  method  of  predicting  solubility  used  in  this  work 
does  not  directly  allow  for  the  simulation  of  a  dynamic  event,  such  as  thermal  rotation,  it  is 
possible  that  this  dynamic  effect  is  indirectly  included  in  this  model  using  the  random  selection 
in  Model  A.  In  NPT-MC  simulations  using  Model  C,  this  dynamic  reorientation  of  the  dipoles  is 
not  incorporated.  Thus  in  Model  C  the  term  <V  exp(-TVkT)>  in  Eq.  (7)  misses  the  dominant 
contributions  of  those  states  of  the  solvent  in  which  the  dipole  moment  vectors  of  the  modifiers 
in  the  neighborhood  of  the  test  particle  are  nearly  parallel  to  each  other  and  antiparallel  to  that  of 
RDX.  This  would  explain  the  success  of  Model  A  in  the  comparison  of  its  predictions  to  those 
of  Model  C. 

(g)  Limitations  of  the  infinite  dilution  approximation: 

The  solubility  of  RDX  in  pure  CO2  at  T=323  K,  P=41.4  MPa  (0.0916  mg  RDX/g  CO2) 
corresponds  to  1  molecule  of  RDX  in  about  55000  molecules  of  CO2.  Such  a  solution  can 
reasonably  be  considered  one  with  infinite  dilution.  The  theoretical  method  used  here  is  invalid 
when  the  solution  is  not  an  infinitely  dilute  solution.  The  calculated  values  of  the  solubility  of 
RDX  in  pure  CO2  indicate  that  for  those  simulations,  the  infinite  dilution  assumption  is  valid.  In 
the  present  calculations,  in  which  we  consider  about  300  molecules  of  solvent  in  the  box,  the 
solubility  is  less  than  1  molecule  of  RDX  in  300  molecules  of  C02  for  SEFs  up  to  180  at  this 
temperature  and  pressure.  Beyond  this  range,  in  the  finite  dilution  case,  one  would  require  a 
simulation  in  the  presence  of  RDX  molecules  and  the  (n+N)  CO2  and  modifier  molecules.  This 
may  alter  the  environment  for  the  test  particle.  The  presence  of  the  RDX  molecules  at  the  low- 
energy-sites  in  the  proximity  of  a  few  modifier  molecules  may  prohibit  sampling  of  such 
favorable  sites  by  the  test  particle,  and  the  test  particle  may  sample  a  larger  number  of  repulsive 
locations  than  in  a  solution  that  does  not  contain  RDX.  These  conditions  would  lead  to  the  lower 
value  of  solubility  than  that  given  by  the  infinite  dilution  approximation. 
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In  view  of  the  foregoing  discussion  it  is  very  likely  that  the  results  for  DMSO  at  the  5% 
and  6%  concentrations  as  reported  in  Table  2  give  an  upper  limit  to  the  solubility,  i.e.,  it  is  likely 
that  at  these  concentrations,  lower  solubilities  than  those  given  in  Table  2  will  result  from 
calculations  using  an  improved  formulation  that  does  not  assume  the  infinite  dilution 
approximation. 

TNT  solubility  in  pure  SF  C02 

The  solubility  of  TNT  in  pure  SF  C02  using  both  monosite  and  multisite  models  of  the 
interactions  between  TNT  and  C02  as  a  function  of  temperature  and  pressure  have  been 
performed,  and  are  reported  with  the  experimental  values  in  Table  5.  The  values  of  solubility 
(both  experimental  and  theoretical)  are  very  large  compared  to  those  of  RDX.  The  basis  for  the 
differences  in  magnitude  of  solubilities  between  the  two  systems  lies  in  the  value  of  p*  used  in 
Eq.  (6).  For  example,  at  T=323  K,.  P=41.4  MPa,  the  ratio  of  the  computed  solubility  of  TNT  and 
RDX  is  333,  whereas  the  corresponding  ratio  of  p*  is  873;  at  T=308  K,  P=13.8  MPa,  these  two 
ratios  are  827  and  1232,  respectively. 


Table  5:  Solubility  of  TNT  in  SCF  CP2 


Temp.  (K) 

Pressure  (MPa) 

Density  (g/cc) 

Solubility  (mg/g) 

Multisite  Monosite  Expt(a)’(b) 

353.0 

27.6 

42.1 

29.2 

13.8 

| 

1.0 

- 

338.0 

27.6 

| 

37.2 

28.5 

BFTTHi» 

13.8 

| 

3.0 

- 

10.4 

0.278 

0.1 

- 

323.0 

41.4 

31.0 

44.2 

27.6 

22.9 

19.4  ™ 

13.8 

0.720 

8.0 

“ 

10.4 

0.443 

1.0 

- 

308.0 

48.3 

1.007 

26.1 

37.6 

41.4 

0.984 

22.4 

21.4 

27.6 

0.931 

15.5 

13.8 

0.823 

9.1 

5.7 

10.4 

0.784 

7.8 

303.0 

48.3 

1.022 

22.2 

25.3 

— — 

41.4 

HHEUiHH 

31.6 

15.5 

mmm 

27.6 

22.2 

13.8 

TFFF 

13.8 

0.867 

10.4 

0.824 

5.9 

6.9 

5.3 

293.0 

37.9 

9.6 

15.3 

WEEtim 

13.8 

hsh 

8.3 

4.3 

BH 

(a) .  Ref.  (28) 

(b) .  Results  of  Ref.  (29)  has  been  further  analyzed  by  Morris  and  are  listed  in  Ref.(28). 
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The  computed  solubility  of  TNT  at  353  K,  27.6  MPa  corresponds  to  ~  2  molecules  of 
TNT  in  the  simulation  box  containing  216  molecules  of  CO2.  This  suggests  that  the  infinite 
dilution  approximation  assumed  in  the  formulated  used  for  these  calculations  is  not  valid  at  this 
temperature  and  pressure  for  the  TNT  sytsem.  For  more  accurate  results  at  this  temperature  and 
pressure,  we  are  required  to  modify  the  model  by  assuming  the  solvent  consists  of  (n+N) 
molecules,  where  n  denotes  a  small  number  of  TNT  molecules  and  N  denotes  the  number  of  CO2 
molecules.  With  this  consideration,  the  interaction  of  the  test  particle  with  these  n+N  molecules 
should  be  computed  to  evaluate  T'  in  Eq.  (7).  By  repeating  the  calculation  for  different  values  of 
n,  we  can  find  n  such  that  the  computed  solubility  corresponds  to  the  ratio  n/N.  The  solubility 
computed  in  this  fashion  would  not  be  based  on  the  infinite  dilution  approximation  and  as  such 
would  be  more  accurate.  For  most  of  the  data  for  TNT  reported  in  Table  5,  however,  the 
solubility  corresponds  to  n«l ;  as  such,  the  infinite  dilution  approximation  is  expected  to  be 
valid  for  all  such  computations. 


Conclusions 

Models  for  use  in  molecular  simulation  of  crystalline  explosives  and  explosives  in  pure 
and  polar-modified  SF  CO2  were  developed.  Simulations  were  performed  to  assess  the  quality  of 
the  models,  with  the  results  indicating  that  a  significant  capacity  for  predictions  using  these 
models  exists.  Dependence  of  solubility  enhancement  factors  due  to  modifier  type,  dipole 
moment,  concentration,  effective  size,  temperature  and  pressure  were  examined.  The  results 
suggest  that  the  degree  of  solubility  enhancement  has  a  cooperative  dependence  on  the  magnitude 
of  the  dipole  moment,  the  number  of  modifier  molecules  surrounding  the  solute  particles,  and  the 
effective  size  and  strength  of  the  modifier-solute  interaction. 


SUMMARY  OF  THE  FINAL  PRODUCT/DELIVERABLE 

The  finished  product  of  PP-695  consists  of  a  set  of  FORTRAN  77  computer  codes  that 
predict  the  solubility  of  a  nitramine  explosive  in  neat  or  polar-modified  CO2.  The  working 
equations  for  the  predictions  are  based  on  classical  statistical  mechanical  theories,  and  the 
models  are  developed  using  the  fundamental  quantum  mechanical  theories.  Computer 
simulations  using  these  codes  provide  a  prediction  of  the  solubility  of  an  ingredient  of  a 
propellant  as  a  function  of  the  conditions  of  a  proposed  experiment,  i.e.  the  temperature, 
pressure,  and  percent  cosolvent  (if  desired).  The  software  is  available  for  electronic  distribution 
through  the  SERDP  Support  office  or  by  contacting  Dr.  Betsy  M.  Rice  (phone:  410-306-1904,  e- 
mail:  betsvr@arl.mil). 
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Figure  Captions 


Fig.  1:  The  solubility  enhancement  factor  (SEF)  and  extraction  enhancement  factor  (EEF)  in 
solvent  for  P=41.4  MPa,  T=  323  K. 

Fig.  2:  The  solubility  enhancement  factor  dependence  on  the  dipole  moment  pb  with  all  other 
parameters  the  same  as  those  used  for  DMSO  simulations  [P=41.4  MPa,  T=323  K,  3.38%]. 

Fig.  3:  The  average  number  of  modifier  neighbors  of  the  RDX  molecule  as  a  function  of  the 
potential,  T  in  kcal/mol,  experienced  by  the  RDX  molecule  as  a  test  particle  [T=323  K,  P=41.4 
MPa,  6%  DMSO]. 

Fig.  4:  The  pair  distribution  function  for  solvent  particles  in  polar-modified  CO2  solvents  at 
P=41.4  MPa,  T=323  K.  Upper  frame  denotes  solvent  with  a  concentration  of  DMSO=3.38%; 
middle  frame  denotes  the  solvent  with  a  CH3CN  concentration  of  3.99%;  and  lower  frame 
denotes  the  solvent  with  a  CH3OH  concentration  of  5.54%. 

Fig.  5:  Energetically-favored  interaction  potential  V(min.)  (see  text)  as  a  function  of  the 
distance  (in  A)  between  the  centers  of  mass  of  modifier  and  RDX  molecules.  DMSO-RDX  and 
CH3CN-RDX  interactions  are  represented  with  circles  and  triangles,  respectively. 

Fig.  6  :  Two-dimensional  illustration  of  three  modifiers  surrounding  RDX  for  the  DMSO/CO2 
system  (left-hand  portion  of  the  figure)  and  CH3CN/CO2  system  (right-hand  portion  of  the 
figure).  Solid  lines  connecting  the  molecules  denote  most  probable  intermolecular  separations 
(determined  through  computer  simulation  in  the  case  of  modifier-modifier  separations  and  using 
standard  arithmetic  mean  rules  for  RDX-modifier  separations). 
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APPENDIX  A 
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